function outputfile = ratheartmakeregmask( studydir );
% ratheartmakeregmask: make a region of interest for registration

% Usage:
%   outputfile = ratheartmakeregmask( studydir );
% Where:
% studydir = the input T2 directory.
% outputfile = the name of the output mask in the postproc dir.
%
% This function is fully automatic. It does not need to be called by the
% user. It is used internally ratheartregister.m
%


% remove slashes at the end
while ( studydir(end) == '/' || studydir(end) == '\' )
    studydir = studydir(1:(end-1));
end

% find postproc dir.
outdir = [ studydir 'postproc\' ];
if ~exist(outdir, 'dir')
  error( 'makeregmask:invalidoutputdir', [ outdir ': invalid postproc dir. Run first ratheartT2fit!' ]);  
end

% read the r2map, which defines the heart mask.
cd(outdir)
d = dir;
for j=3:size(d,1)
    x=d(j).name;
    if size(char(regexp(x,'R1map', 'match')),2)>=1 %read in the fixed longest IR image
        info = dicominfo( fullfile(outdir, 'R1map.dcm') );
    elseif size(char(regexp(x,'R2map', 'match')),2)>=1 %read in the fixed longest IR image
        info = dicominfo( fullfile(outdir, 'R2map.dcm') );
    end   
end
img = double(dicomread(info));

for dim = 1:numel(size(img))
  % find first and last row that is nonzero, for each dimension
  rows = max(abs(img),[],dim) >0;
  boxl(dim) = find( rows, 1);
  boxu(dim) = numel(rows) + 1 - find( rows(end:-1:1), 1);
end
% enlarge box a bit.
boxsize = boxu - boxl +1;
boxl = round(boxl - 0.5*boxsize);
boxu = round(boxu + 0.5*boxsize);
% make sure it's not bigger than the original image
boxl  = max(  boxl , [1 1] );
boxu  = min(  boxu , size(img) );
boxsize = boxu - boxl + 1;
% generate the mask
imagemask = zeros( size(img) );
% confuses x/y dimensions:
% imagemask( boxl(1):boxu(1), boxl(2):boxu(2) ) = ones(boxsize);
% very confusing, but in dicom: x=horiz,y=vert, matlabmatrix x=vert,y=horiz
imagemask( boxl(2):boxu(2), boxl(1):boxu(1) ) = ones(boxsize)'; %


% debug:
%figure; imshow(img);
%figure; imshow(imagemask);

% save the mask
outmask = fullfile(outdir, 'regmask.dcm' );
dicomwrite( int16(imagemask), outmask, info );

outputfile = outmask;
